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Abstract. This article describes a comparison of two calculations of the 
merger of a binary neutron star (NS) system which is initially within the 
tidal instability as described by Rasio & Shapiro [0, ||. The same initial 
data is used with one simulation involving a purely Newtonian evolution 
of the Euler equations for compressible fluids and another similarly evolved 
except there is an inclusion of a gravitational radiation reaction (GRR) term 
at the 2.5 Post-Newtonian (PN) order as prescribed by Blanchet, Damour, 
& Schafer The inclusion of GRR is to allow an approximation of the 
full relativistic effect which forces the inspiral of binary systems. The initial 
data is identical at the start of each evolution and only the inclusion of 
the 2.5PN term differs in the evolutions. We chose two co-rotating, r=2 
polytropic stars with masses IAMq, radii i?=K=9.56km and central densities 
of 2.5 X 10^^ g/cm^. The initial binary separation is 2.9-R*, inside the region 
where the dynamical tidal instability has been shown to cause the merger 
of the two stars. Comparisons between measured quantities will indicate a 
substantial difference between the two evolutions. It will be shown that the 
effect of GRR on the dynamics of the merger is signiflcant. 



1. Introduction 

We study the evolution of binary NS systems as it represents one of the most interesting 
astrophysical systems of the known universe. Since the discovery of the flrst binary pulsar 
system in 1974 by Hulse & Taylor 0] and the subsequent measurement of the orbital 
period decrease as predicted by general relativity it is has been thought that the energy 
loss from these systems, in the form of gravitational waves (GWs), could be measured. 
Such measurements are expected to begin early in the next century when several GW 
observatories or detectors, such as the US-based LIGO experiment 0, are to come on- 
line in their advanced stages. These systems also represent interesting laboratories for 
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nuclear physics. Most recently, the site of 7-ray bursts has been suggested as binary NS 
mergers^, ^. See Calder, Swesty & Wang 1988 for discussion flOf . 



2. Method of simulation 



2.1. Numerical hydrodynamics and self gravity 



We assume the neutron star can be described by a compressible fluid which satisfies the 
Euler equations, plus the solution to Poisson's equation for self-gravity and the 2.5PN 
GRR potential ^ which enters into the momentum equation (see [§, |l^ for a discussion 



of the the Euler equations with self-gravity). The method is similar to the one employed 
in the ZEUS codes [|ll[]. Usage of the 2.5PN GRR potential is consistent with for 



hydrodynamic evolution up to 2.5PN order. 

2.2. 2. 5PN gravitational radiation reaction potential ^ & boundary conditions 

Solving for the 2.5PN GRR potential ^ at any time during the evolution requires a 
solution of a Poisson equation of the form 

V'i? = A-kGD'^x,^ (1) 

wheYeD^^=d^i^'' /dt^ is third time-derivative of the symmetric, trace- free mass quadrupole 
tensor, and the 2.5PN GRR potential is defined as 

We utilize a full multi-grid W-cycle algorithm to calculate the solution of (1). To define 
the boundary values of R we utilize a Riemann sum decomposition similar to what is 
described in for the Newtonian potential. The solution to (1) will be of the form, 

R=-G f'^ dr' (3) 



where f is the distance from the mass element to the boundary zone (see [jTO[ for dis- 
cussion). We have found that R=0 boundary conditions are sufficient. 

3. Simulations 

3.1. Technical concerns: Evolution time scales 

The high densities present in neutron stars implies the adiabatic sound speed Cs will 
be very large. The Courant-Friedrichs-Lewy (CFL) condition defines the minimum time 
step and since it is dependent on the inverse of Cg one expects evolutions on physical 
time scales (milliseconds) to require many thousands of numerical time steps. Explicitly, 
the CFL condition may be written as 

At = a (Ax/Cs) (4) 

where a is the Courant number. The time step may also be reduced if high velocities 
are present or if there is dissipation A suite of code tests, spatial and temporal 



resolution tests, and questions of stability have been addressed [|T0] . 
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3.2. Initial data and atmosphere 



We construct initial data for a close binary system because of the computational imprac- 
ticality of evolving two stars from infinity. Initially, two non-rotating, r=2 polytropic 
stars with masses I.AMq, radii -R*=9.56 km and central densities of 2.5 x 10^^ g/cm^ 
are placed on a grid with an initial binary separation of 2.9R^: between the star centers. 
The equation of state we chose is a polytropic one which involves a relationship between 
pressure P and density p, only. It should be noted that the equation of state for the 
evolutions is not isentropic but is an ideal gas equation of state defined in terms of 

P = {T -1)E = {T -l)ep (5) 

where e is the specific internal energy. We prescribe a low density atmosphere, p=1.0 x 10^ 
g/ cm'^, with a high internal energy. This prevents the atmosphere from falling onto the 
star surface and shocking over short time scales. We have found other treatments of 



the atmosphere, such as the one employed by Ruffert et al |T2[, to be problematic and 
erroneous. All the fluid elements in the stars and a nearby region of atmosphere are then 
set to co-rotate about the center of the grid with a frequency defined by Kepler's third 
law for two point particles. The axis of rotation is the z-axis. 



3.3. Results 

For the two simulations we consider, the grid size is 129^ with a (ix=0.89 km allowing for 
a resolution of ~21 zones across each star. The Newtonian evolution, simulation N2.9, 
is evolved for 3776 time steps with a total physical time of 5.40ms. Tidal distortion of 
the stars begins immediately. By one-eighth of an orbit the surfaces of the stars are 
touching. After one-half an orbit only the densest regions of the stars are not in contact. 
However, it is after this time during which the cores of the individual stars begin to 
fall inward toward the system center of mass that we find mass shedding and angular 
momentum transfer from the cores to the lower density regions. As these low density 
regions gain angular momentum spiral arms form carrying this material away from the 
coalescing object and off the evolution grid. By seven-eights of an orbit the densest parts 
of the stars are touching and the spiral arms are pronounced. Soon after the peak in the 
gravitational wave luminosity (GWL) occurs the magnitude of the GWL and the GWs 
h-\. and hy, falls-off dramatically, as depicted by /ix in figure 1. A differentially rotating 
object remains by one and a half orbits with the merged core rotating faster than the 
lower density disk surrounding it. The core is non-axisymmetric, as would be indicated 
by a non-zero waveform h-^ (figure 1). 

The PN2.9 simulation where GRR at 2.5PN order is included is evolved for 3694 
time steps with a total physical time of 5.26ms. The results are qualitatively similar to 
the Newtonian simulation except that the merger occurs substantially earlier. 

The simulations have some substantial differences in the gravitational waveforms, 
the time scales to merger, and the equilibrium conserved quantities. A phase difference 
develops in h-^ after 0.8ms due to the inclusion of the GRR which accelerates the merger. 
By 1.5ms, a 90 degree phase difference can be seen in /ix! the merger has occured a 
quarter orbit earlier. The mass loss that occurs is the means by which the pre-merger 
stars shed angular momentum to form the central object. Since GRR reduces the total 
angular momentum of the system throughout the evolution in PN2.9 we see a higher total 
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mass of the post-merger object or less mass shedding than in the N2.9 calculation. Both 
mass and angular momentum are less in the post-merger N2.9 object after the merger 
(figure 2). The inclusion of GRR alters the evolution of the system and the structure of 
the post-merger object. Exclusion of GRR clearly ignores an effect which significantly 
changes the evolution of a binary NS system. This effect is not overwhelmed by the 
Newtonian dynamical tidal instability of ||l|, 0. 

At 3.485ms into the evolution and nearly 1.8ms past the merger the core region can 
be defined as having a torus-like mass distribution (figure 3). Outward from the center 
to about 3km the mass density rises along the x- and y-axes then falls off rapidly. As 
seen in figure 4, most of the mass in the system is contained within 27km of the system 
center of mass. However, a significant fraction of the kinetic energy of the system lies 
outside this region (figure 5). The merged object is differentially rotating (figure 6) with 
the merged core rapidly rotating inside low density disk. The post-merger object appears 
to be nearly axisymmetric, and more so than in N2.9. 

In conclusion, the effect of GRR on the evolution of the binary NS merger is signifi- 
cant. There is a substantial difference between the Newtonian and Newtonian-|-2.5PNRR 
evolutions. The final merged objects are different in structure and total angular momen- 
tum. The differing pre-merger inspiral rate produces a phase difference in the waveforms 
and the total radiated energy through GWs is somewhat different. Thus, it can be 
stated that the effect of the 2.5PNRR potential is significant in the late portion of the 
merger and in forming the final coalesced object. This further motivates the need for 
fully relativistic simulations. 
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Figure 1 Figure 2 
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Figure 1. Gravitational waveform /ix is shown for 5.26ms as seen by an observer along the 
z-axis. 

Figure 2. The mass for a similar time at the four density thresholds. 
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Figure 3. Mass density as a function of radius. 

Figure 4. Total mass enclosed as a function of radius as a solid line while the dotted lines are 
factors of the Schwarzschild radius Rs- 

Figure 5. Kinetic energy as a function of radius. 

Figure 6. The norm of the fluid velocity as a function of the radius. 
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